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All means (even continuous) sanctify the discrete end. Doron Zeilberger 2 

Abstract: The n th Birkhoff polytope is the set of all doubly stochastic n x n matrices, that is, those matrices with 
nonnegative real coefficients in which every row and column sums to one. A wide open problem concerns the volumes 
of these polytopes, which have been known for n < 8. We present a new, complex-analytic way to compute the 
Ehrhart polynomial of the Birkhoff polytope, that is, the function counting the integer points in the dilated polytope. 
One reason to be interested in this counting function is that the leading term of the Ehrhart polynomial is — up to a 
trivial factor — the volume of the polytope. We implemented our methods in the form of a computer program, which 
yielded the Ehrhart polynomial (and hence the volume) of the ninth and the volume of the tenth Birkhoff polytope. 

1 Introduction 

One of the most intriguing objects of combinatorial geometry is the n th Birkhoff polytope 
' ( xn ■■■ x ln 

I (Z jj n2 ■ x ■ ^ 

3 ~ ' X^fe x jk = 1 f° r all 1 < j < n 



Y^j Xjk = 1 for all 1 < k < n 

■ \ - ■"■ - ' v. 

, \ X n i 



often described as the set of all n x n doubly stochastic matrices. B n is a convex polytope with 
integer vertices. It possesses fascinating combinatorial properties [4, 5, 6, 8, 23] and relates to 
many mathematical areas [10, 14]. A long-standing open problem is the determination of the 
relative volume of B n , which had been known only up to n = 8 [7, 17]. In this paper, we propose a 
new method of calculating this volume and use it to compute volSg and vol£>io- 

One of the recent attempts to compute vol B n relies on the theory of counting functions for the 
integer points in polytopes. Ehrhart proved [11] that for a polytope V C W 1 with integral vertices, 
the number 

L v (t) := # (tvnz d ^ 



is a polynomial in the positive integer variable t. He showed various other properties for this 
counting function (in fact, in the more general setting of V having rational vertices), of which we 
mention three here: 



• The degree of L-p is the dimension of V . 

• The leading term of L-p is the relative volume of V, normalized with respect to the sublattice 
of 7L d on the affine subspace spanned by V. 

1 Appeared in Discrete & Computational Geometry 30, no. 4 (2003), 623-637. 
2 E-mail to the first author on Wed, 14 Aug 2002 16:29:35 -0400 (EDT) 
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Since L-p is a polynomial, we can evaluate it at nonpositive integers. These evaluations yield 

L P (0) = x(V), (1) 
L P (-t) = (-lf^L P o(t) . (2) 



Here x('P) denotes the Euler characteristic, V° the relative interior of V. The reciprocity law (2) 
was in its full generality proved by Macdonald [15]. 

The application of this theory to the Birkhoff polytope B n incorporates the nice interpretation of 
the number of integral points in tB n as the number of semi-magic squares, namely, square matrices 
whose nonnegative integral coefficients sum up to the same integer t along each row and column. 

We will denote the Ehrhart polynomial of B n by 

H n {t) :=L Bn (t) . 

It is not hard to see that dimB n = (n — l) 2 , hence H n is a polynomial in t of degree (n — l) 2 . The 
first two of these polynomials are trivial: 

#i(t) = l, H 2 (t)=t + 1, 

the first nontrivial case was computed by MacMahon [16] as 

^-(TMT)- < 3 > 

The structural properties of H n were first studied in [12, 18, 20]. It is a nice exercise to deduce 
from (2) that 

H n (-n-t) = (-l) n - 1 H n (t) (4) 

and 

H n (-1) = H n (-2) = ■■■ = H n (-n + 1) = . 

This allows the following strategy of computing H n , and therefore, the volume of B n : compute 
the first ( U 2 1 ) values of H n , use the above symmetry and trivial values of H n , and calculate the 
polynomial H n by interpolation. In fact, as far as we are aware of, the volume of Bs was computed 
using essentially this method, combined with some nice computational tricks [7, 17]. 

We propose a new, completely different approach of computing H n (and hence vol£>„). It is based 
on an analytic method by the first author of computing the Ehrhart polynomial of a polytope [2]. 
We will introduce the application of this method to the Birkhoff polytope in the following section. 

Some recent refreshing approaches — of a more algebraic-geometric/topological flavor — to the prob- 
lem of computing vol0 n and H n can be found in [1, 9, 22]. 
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2 An integral counting integers 



A convex polytope P C l rf is an intersection of halfspaces. This allows the compact description 

V = |x G R d : Ax < b j , 

for some (m x (i)-matrix A and m-dimensional vector b. Here the inequality is understood compo- 
nentwise. In fact, we may convert these inequalities into equalities by introducing 'slack variables.' 
If V has rational vertices (those polytopes are called rational), we can choose A and b in such a 
way that all their entries are integers. In summary, we may assume that a convex rational polytope 
V is given by 



V 



-{ 



x G 



v >0 



: Ax 



(5) 



where A G M mX( ^(Z) and b G Z m . (If we are interested in counting the integer points in V, we may 
assume that V is in the nonnegative orthant, i.e., the points in V have nonnegative coordinates, as 
translation by an integer vector does not change the lattice-point count.) The following straight- 
forward theorem can be found in [2]. We use the standard multivariate notation v w := v™ 1 ■ ■ ■ v™ n . 



Theorem 1 [2, Theorem 8] Suppose the convex rational polytope V is given by (5), and denote the 
columns of A by Ci , . . . , . Then 



Lr(t) 



1 



-t&i-i 



{2my- j\ Zl \ =ei J\z m \=e r , 
Here < e± , . . . , e m < 1 are distinct real numbers. 



(1 



)---(i-z°«; 



dz 



It should be mentioned that L-p is in general not a polynomial if the vertices of V are not integral, 
but a quasipolynomial, that is, an expression of the form 

c d (t)t d + • • • + a(t)t + co(t) , 

where Co, . . . , are periodic functions in t. (See, for example, [19, Section 4.4] for more information 
about quasipolynomials.) Theorem 1 applies to these slightly more general counting functions; 
however, in this article we will only deal with polytopes with integer vertices, for which L-p is a 
polynomial. 

We also note here that Theorem 1 can be used to quickly compute by hand formulas for certain 
classes of polytopes (see, for example, [2]). In this project, we take a slightly different approach 
and use this theorem to efficiently derive formulas with the help of a computer. 

We can view the Birkhoff polytope B n as given in the form of (5), where 

( 1 ••• 1 \ 
1 ... i 



A = 



1 



1 I 



3 



is a (2n x n 2 )-matrix and b = (1, . . . , 1) € l? n . Hence Theorem 1 gives for this special case 

hm) = 1 f ... f {zi---z 2n y t - 1 dz 

(2TTl) 2n J J (1 - ZlZn+l)(l - Z 1 Z n+2 ) ■ ■ ■ (1 - Z n Z 2n ) 

Here it is understood that each integral is over a circle with radius < 1 centered at 0; all appearing 
radii should be different. We can separate, say, the last n variables and obtain 

Hn{t) = J^— f ■■■ [(Zl--- Zn)- 1 - 1 f Z - dz] dz n --- dz X . 

(2iri) n J J \2iri J (1 - z x z) ■ ■ • (1 - z n z) J 

We may choose the radius of the integration circle of the innermost integral to be smaller than the 
radii of the other integration paths. Then this innermost integral is easy to compute: It is equal 
to the residue at of 

1 



z t + 1 (l-z 1 z)---(l-z n z) 

and, by the residue theorem, equal to the negative of the sum of the residues at z-f 1 , . . . , z^ 1 . (Note 
that here we use the fact that t > 0.) The residues at these simple poles are easily computed: the 
one, say, at 1/zi can be calculated as 

1 \ 1 z ~ z~ 1 

lim ( z —r—r, ; ; r = lim 



z->l/zi V Z\) Z t+1 (l — Z\Z) ■••(]. — Z n z) z->l/zi 1 — Z\Z z t+1 (l — z 2 z) ■••(! — z n z) 
1 „t+l jt+n-l 



z\ (1 - z 2 /zi) • • • (1 - z n /zi) {zi - z 2 ) ■ ■ ■ Oi - z n ) 

This yields the starting point for our computations. 
Theorem 2 For any distinct < e±, . . . , e n < 1, 

1 r r ( n z* +n_1 \ n 

(2iri) n J ]zil=ei J ]Znl=en Uj^k( Zk - z ^ J 

Remark. It can be proved from the form of the integrand that H n is indeed a polynomial in t: To 
compute the integral, one has to execute a (huge) number of limit computations, which yield "at 
worst" powers of t (as a consequence of L'Hospital's Rule). In fact, one can make this property 
more apparent by noticing that the expression in parenthesis is actually a polynomial; namely 

(27™) J\*i\=*i J \*n\=e n \ mi+ ±? mn=t J 

where the sum is over all ordered partitions of t. This formula can also be proved "more directly" 
combinatorially. 3 

3 The authors thank Sinai Robins and Frank Sottile for their help in the proof of this equivalence and its combi- 
natorial interpretation. 
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3 Small n do not require a computer 



We will now illustrate the computation of H n (and hence vol B n ) by means of Theorem 2 for n = 3 
and 4. These calculations "by hand" give an idea what computational tricks one might use in 
tackling larger n with the aid of a computer. 

By the theorem, 

If ( z t+2 z t+2 z t+2 \ 3 

H 3 (t) = 777-73 / (zizazs)"*" 1 t h v + t h 7 + 7 y, 7 dz. 

(2-Kiy J \(z l - z 2 ){z 1 - z 3 ) (z 2 - z 1 )(z 2 - z 3 ) (z 3 - Zl){Z 3 - Zi) ) 

We have to order the radii of the integration paths for each variable; we choose < e 3 < e 2 < e± < 1. 
We heavily use this fact after multiplying out the cubic: integrating, for example, the term 

~-t-i _-t-i ^2t+5 



(z 3 - z 2 ) 3 (z 3 - Zl] 



3 



with respect to z 3 gives 0, as this function is analytic at the ^-origin and \z\ \ , \z 2 \ > e 3 . After 
exploiting this observation for all the terms stemming from the cubic, the only integrals surviving 
are 

2t+5 -t-l -t-l 

1 " 2 Z3 dz 



(2vn)3 J ( Zl - Z2 f( Zl - z rf 
and 

3 f z\+ 3 z^- 1 



I 



dz . 



(2m) 3 J (zi - z 2 ) 3 (z 1 - z 3 ) 2 (z 2 - z 3 ) 
The first integral factors and yields, again by residue calculus, 

1 f zf^z^z^' 1 1 f 2t+5 ff z-t- 1 

(2m) 3 J ( Zl - z 2 f{ Zl - z 3 ) 3 dZ - (2m) 3 J ^ \J ( Zl - z) 3 dZ ) dZl 

For the second integral, it is most efficient to integrate with respect to z 2 first: 

v t+3^^-t-i q r j+3 z -t 

dz 3 dz\ 



3 f zjr^z^ 3_ f gTj 

(2m) 3 J ( Zl - Z2 f( Zl - Z3 y( Z2 - Z3 ) (2mfj ( Zl - z 3 f 

= ~ J 4 +3 ^(-t)(-t ~ !)(-< - 2)H " 3) z?- 4 dz, = -3( 



t + 3 
4 



Adding up the last two lines gives finally 



which is equal to (3). To obtain the volume of B 3 , the leading term of H 3 has to be multiplied by 
the relative volume of the fundamental domain of the sublattice of Z 9 in the affine space spanned 
by B 3 . This volume is 9; hence 

vol #3 = - . 



5 



In general, it is not hard to prove (see, for example, the appendix of [7]) that the relative volume 
of the fundamental domain of the sublattice of Z n in the affine space spanned by B n is n n_1 . 

The number of integrals we have to evaluate to compute H 4 is only slightly higher. By Theorem 2, 

4 

I dz . 



i+3 



= ttm I I ! I i^mr 1 - 1 {jZvr4 

(2vr*) 4 J ]zi]=ei J\ Z2 \ =e2 J\ Z3 \ =e3 J\ Z4 \ =e4 UjM Zk - z i>. 

Again we have a choice of ordering the radii; we use < £4 < £3 < e 2 < e\ < 1. After multiplying 
out the quartic, we have to calculate five integrals; their evaluation — again straightforward by means 
of the residue theorem — is as follows. As before, we can 'save' computation effort by choosing a 
particular order with which we integrate and by factoring an integral if possible. 

(27Ti)4 J ( Zl - Z 2 y( Zl - Z,)\ Zl - Z,Y ^ - (2vTz) 4 J ^ [J ( Zj - Z f ** ) ^ 

v 2t+8„2 1 



/ 



(27TZ) 4 J (zi - z 2 ) 4 {zi - z 3 ) 3 (zi - z 4 ) 3 (z 2 - z 3 )(z 2 ~ Z4) 



dz 



2M-8? / r „-t-i \ 2 



1 1 ' 1 dz I dz\ dz 2 



(2vri) 4 J ( Zl - z 2 f \J ( Zl - z) 3 (z 2 - z) 

z x z 2 (Jt + 2\z l +2{t+1) z 1 +2 ^i _*? ) dz 



(2m) 2 J (zi-z 2 ) 7 \ \ 2 Jz 1 -z 2 " (z 1 -z 2 ) 2 (z 1 - z 2 ) 3 (zi - z 2 ) 3 

"t + 2\/t + 5\ n . _/t + 6\ /* + 7\ /2i + 8\ 



I 



r) + 8, (+1 ,er) +8 (- 7 )-e 

^t+S^-t-l^y-t-l 

dz 



9 J 



(27ri) 4 J (zi - z 2 ) 3 (zi - z 3 ) 4 0i - z 4 ) 3 (z 2 - z 3 )(z 3 - z 4 ) 



_2t+8 -t+1 -t-1 

4,1 6r) & A 

dz 



{2m) 3 J (zi - z 2 ) 3 (z 1 - z A ) 7 (z 2 - z 4 ) 

z x z, ((t + 2\ Zl + + - + -J1 -)dz 



(2ni) 2 J (z\ — Z4) 7 \\ 2 )z\ — z 4 (z\ — z 4 ) 2 (z\ — z 4 ) 

<cr)cr)+<-»cr)+cr 

(■ f yt+5 t+5 _-t~l y-t-1 



(2vTz) 4 J ( Zl - Z 2 )\z X - Z 3 ) 2 { Zl - Z 4 ) 2 (z 2 - Z 3 ) 2 {z 2 - Z 4 ) 2 

t+5j+5 / r y-t-\ \ 2 



6 f z^°zV° f z 



(2mrJ (z x - Z2 y \) ( Zl - z ) 2 ( Z2 - z ) 2 



dz dz\ dz 



2 



a r yt+5 t+5 / _-2t-4 _-2t-3 _-2t-2 



(2vri) 2 J ( Zl - z 2 y V ' (zi ~ ziY (zi- z 2 f ( Zl - z 2 f 



dz 



6 



12 f 7 t+b 7 2 7 2 

z l Z 2 Z 3 Z 4 ^ z 



(2m) A J (zi - z 2 ) 3 (zi - z 3 ) 3 (zi - z A ) 2 (z 2 - z 3 ) 2 (z 2 - z^)(z 3 - z 4 ) 

19 r r t+5„2 r -t+l 

' Zl Z2Z * dz 



{2irif J ( Zl - Z2 f( Zl - ZA f( Z2 - ZA )Z 

12 f ^+V +1 / 1 z 4 zj_ 

(2m) 2 J ( Zl - z 4 ) 5 \(zi - z,Y + % - z 4 ) 4 (zi- z,f 1 ' Z 

-< 5 )-< ( r)-<T 

Adding them all up gives 

„ , , 11 „ 11 „ 19 7 2 fi 1109 . 43 4 35117 , 379 , 65 

F 4 (t) = t 9 H t -\ £' H — t 6 H t 5 H i H t 3 H t 2 H t + 1 

4K ' 11340 630 135 3 540 10 5670 63 18 

and hence 

^ 11 176 
vol£ 4 = 4 3 



11340 2835 



4 Larger n do 



As we have seen in the examples, after multiplying out the integrand of Theorem 2, many of the 
terms do not contribute to the integral. The following proposition will provide us with a general 
statement to that effect. 

For a rational function / in n variables Zj we use the notation d r (f) for the degree of / in the 
variables z\ , . . . , z r . 



Proposition 3 Suppose p±,...,p n are integers, qjk are nonnegative integers (1 < j < k < n), 
1 > ei > • • • > e n > 0, and 



Ui< 



P3 



f( z , z ) - t±l^izi 

If 1 < r < n and d r (f) < —r then 

I ■■■/ /(z)dz = 0. 



Proof. We need only show that 



/ •••/ /(z)dz = 0, (*) 

J\zi\=ei J\z r \=e r 



for then we can integrate over all n variables by first integrating over z\...z r . 
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If r = 1 then /, considered as a function of z±, is analytic outside the circle \z±\ = e± and has zero 
residue at infinity since its degree is less than —1, so (*) follows. 

We continue by induction, so suppose r > 1 and (*) is true for smaller r. Suppose that d r (f ) < — r; 
we may assume d r -i(/) > —{ r — !)• Note that d r {f) > d r ^i(f) + d where d is the degree of 
/ in the single variable z r (the discrepancy is the sum of the exponents qj r for j < r). Hence 
d < d r (f) — d r -i(f) < — 1. We consider / as a function of z r and apply the residue theorem to 
the region outside the circle \z r \ = e r . As above / has zero residue at infinity, so we only need to 
consider the residues at the poles Zj for j < r. Evaluating these residues converts the integral of / 
into a (possibly huge) linear combination of integrals of functions of the same form as /, but in the 
n — 1 variables zi, . . . , z r -±, z r+ ±, . . . ,z n . If g is any one of these functions then we easily calculate 
d r -i{g) = d r (f) + 1 < — (r — 1), and, by the induction hypothesis, the integral of g is zero. □ 

From this proposition we obtain the starting point for our 'algorithm.' 

Corollary 4 For 1 > e x > • • • > e„ > and t>0, H n {t) = 

i r r / \ n f v t+n-i \ mk 

(dW -/ f*"'^-'- 1 E* L " m )n n V;) dz ' 

^ m > J\zi\=€i JW\=e n mi+ ... +mn=n \m 1 ,...,m n j k=i \ll jftk (z k Zj)j 

where ^2* denotes that we only sum over those n-tuples of non-negative integers satisfying m\ + 
• • • + m n = n and mi + • • • + m r > r if 1 < r < n. 



Remark. The condition on mi, ... , m n can be visualized through lattice paths from (0, 0) to (n, n) 
using the steps (1, mi), (1, 7712), . . . , (1, m„). The condition means that these paths stay strictly 
above the diagonal (except at the start and end). 4 

Proof. By Theorem 2, 

t+n-l 



- <sW ■ ■ ■ fa ■ ■ ■ -»-'-' X L " , J n (tS^)) ~ dz • 



- dz 

t+n-l 



We select a partition mi + • • • + m n = n and rewrite the corresponding integrand in the language 
of Proposition 3: 

Pj = (n — l)mj + t(mj — 1) — 1 , qjk = m,j + m^ . 
Now suppose 1 < r < n. The degree of the denominator in zi, . . . , z r is 

r n r n r r 

(mj + rrifc) = (n — 1) mj + r m^ = (n — 1) + nr — r 

j=l fc=j+l j=l j=r+l j=l j=l 



= (n — 1) mj + nr — r ^^(m^ — 1) — 



,2 



4 The authors thank Lou Billera for pointing out this lattice-path interpretation. 
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We subtract this from Y^j=iPj t° get 



r 



d r (f) = {t + r) J2(mj - 1) - r(n - r) 



— r 



Here t is non-negative and r(n — r) is positive, so Proposition 3 implies that the integral is zero 



Theoretically, Corollary 4 tells us what we have to do to compute H n . Thus the integrals we 
computed in our calculations of H% and H4 are exactly the non-zero integrals according to this 
result. For practical purposes, however, the statement is almost worthless for larger values of n. 
The first problem is that the number of terms in the sum equals the (n — l) th Catalan number 



which grows exponentially with n. Another slippery point is the evaluation of each integral. As 
we have seen in the examples, and as can be easily seen for the general case, we can compute 
each integral step by step one variable at a time. However, this means at each step we convert 
a rational function into a sum of rational functions (of one variable less) by means of the residue 
theorem. Again, this means that the number of (single-variable) integrals we have to compute grows 
immensely as n increases. In fact, if we just 'feed' the statement of Corollary 4 into a computer and 
tell it to integrate each summand, say, starting with zi, then Z2, and so on, the computation time 
explodes once one tries n = 7 or 8. We feel that computationally this is as involved as calculating 
a sufficient number of values of H n and then interpolating this polynomial. However, complex 
analysis allows us some shortcuts which turn out to speed up the computation by a huge factor and 
which make Corollary 4 valuable, even from a computational perspective. These 'tricks' all showed 
up already in the examples and include 

1. realizing when a function is analytic at the Zfc-origin, 

2. trying to choose the most efficient order of variables to integrate (based on estimating how 
many terms will be generated by the residue calculations, for each available variable), and 

3. factoring the integral if some of the variables appear in a symmetric fashion. 
Finally, if we are only interested in the volume of B n , we may also be 

4. suppressing a particular integral if it does not contribute to the leading term of H n . 

It is worth noting that each of these computational 'speed-ups' decreases the total computation 
time substantially. By applying them to Corollary 4, we implemented a C++ program for the 
specific functions we have to integrate to compute H n . We were able to verify all previously known 
polynomials (n < 8) and to compute vol£>g, vol Bio, an d Hg. The results of our calculations, 
including the polynomials H n for n < 9, and the source code for our program are available at 
www.math.binghamton.edu/dennis/Birkhoff. The following table gives the volumes together 
with the respective computing time (on a 1GHz PC running under Linux). Note that computing the 



unless J2 T j=i( m j — 1) > 0. 



□ 
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full polynomial H n takes longer, as we cannot make use of shortcut #4. In fact, the computation of 
Hg took about 325 days of computer time, although the elapsed time was only about two weeks since 
we distributed the calculations among a number of machines (8 - 40, depending on availability). 



1 b 


vol B 


1 11X16 


1 


1 


< .01 sec 


2 


2 


< .01 sec 


3 


9 


< .01 sec 


4 


176 

6000 


< .01 sec 


5 


23590375 
167382319104 


< .01 sec 


6 


9700106723 
1319281996032-10 6 


.18 sec 


7 


77436678274508929033 
137302963682235238399868928-10 8 


15 sec 


8 


5562533838576105333259507434329 
12589036260095477950081480942693339803308928-10 10 


54 min 


9 


559498129702796022246895686372766052475496691 


317 hr 


92692623409952636498965146712806984296051951329202419606108477153345536-10 14 



After this article was submitted for publication, we used the idle time on our departmental machines 
to compute the volume of B\q. After a computation time of 6160 days, or almost 17 years (again 
scaled to a 1GHz processor), we obtained 

vol #10 = 

727291284016786420977508457990121862548823260052557333386607889 
828160860106766855125676318796872729344622463533089422677980721388055739956270293750883504892820848640000000 ' 

We hope that the reader will take this immense computing time as a challenge to improve our 
algorithms, and work towards n =11. 



5 An outlook towards transportation polytopes 

The Birkhoff polytopes are special cases of transportation polytopes, which are defined below. The 
study of this class of polytopes, which are naturally at least as fascinating as the Birkhoff polytopes, 
was motivated by problems in linear programming; for combinatorial properties see, for example, 
[10, 13]. The goal of this section is to show how our methods can be applied in this more general 
setting. 

Fix positive real numbers a±, . . . , a m , bi,...,b n such that a\ + ■ ■ ■ + a m = &! + ••• + b n . Let 
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a = (ai, . . . ,a m ),b = (61, . . . ,&„). 

(Zll • ■ ■ Xi n \ 
: . el . x jfc >0, Eja . jfc = 6fcfaraU1 < fc < n 

X m l . . . x mn J 

is the set of solutions to the transportation problem with parameters a, b. It is a convex polytope 
of dimension (m — l)(n — 1) in W nn — hence we refer to T a ^ as a transportation polytopes. Note 
that B n = Ti i where 1 = (1, . . . , 1) G Z n . We will be exclusively interested in the case when 
ai, . . . , a m , 61, . . . , b n are integers. One reason for this is that the vertices of 7^ b are then integral. 
Let 

T(a, b) = T(ai, . . . , a m , h, . . . , b n ) = # (T a , b D 

denote the number of integer points in the transportation polytope with integral parameters a, b. 
Geometrically, each of these parameters determines the position of a hyperplane bounding the poly- 
tope T^ h- It is well known that T(a, b) is a piecewise-defined polynomial in a\, . . . , a m , b\, . . . , b n . 
The regions in which T(a, b) is a polynomial depends on the normal cones of the polytopes involved; 
for details, we refer to [21]. 

Moreover, one can derive results for the counting function T(a, b) which are 'higher-dimensional' 
extensions of (4) , by application of the main theorem in [3] , which in turn is a generalization of the 
Ehrhart-Macdonald reciprocity theorem (2). To state this theorem, we need to define the following 
integer-point counting function. Suppose V is a rational polytope given in the form (5); denote by 
"P(t) = {x € M. d : Ax = t} a polytope which we obtain from V = V(b) by translating (some of) 
its bounding hyperplanes. (Classical Ehrhart dilation is the special case t = th.) Let 



£ P (t) = #(p(t)nzj . 



[3, Theorem 4] states that C-p and C-p° are piecewise-defined (multivariable) quasipolynomials 
satisfying 

£ v (-t) = (-l) dimV £ po (t) . 

As easily as (4) follows from (2), this reciprocity theorem yields a symmetry result for the trans- 
portation counting function. Let Id denote the d-dimensional vector all of whose entries are one. 
Then 

T(-a - nl m , -b - ml n ) = (-1)"— 1 T(a, b) . 

We finally turn to the problem of writing T(a, b) in form of an integral. As with B n , we can view 
T^b as given in the form of (5) and apply the philosophy of Theorem 1 to obtain 

T(a, b) = — — / • • • / — ^r 21 77 1 s — 1 dz dw . 

i<k< n J 

Again the integral with respect to each (one-dimensional) variable is over a circle with radius < 1 
centered at 0, and all appearing radii are different. As in the Birkhoff case, we can separate, say, 
the w-variables and obtain 

T(a, b) = — — / • • • / z7 ai - x ■ ■ ■ zZ?™- 1 TT / FFm ^ \ dw k dz . 

v ' ; (27ri)™+«y j 1 m nr=i(i-^») 
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And as with the Birkhoff polytope, the innermost integrals are easy to compute by means of the 
residue theorem. This yields the following statement which generalizes Theorem 2. 

Theorem 5 For any distinct <C e^, . . . , e m 

< 1, 

-iff n m z b k +m-l 
T(a,b) = — - i— / ••• / zr ai_1 • • • O™" 1 TT V r dz . 

Remark. As with Theorem 2, it can be proved from the form of the integrand that T(a, b) is indeed 
a piecewise-defined polynomial in a±, . . . , a m , b\, . . . , b n . 

Theoretically we could now use this theorem to produce formulas for T(a, b) just as we did for 
H n (t). There is one major difference: T(a, b) is only a piecewise-defined polynomial. In fact, we can 
see this from the form of the integral in Theorem 5: whether we will get a nonzero contribution at a 
certain step in the computation depends heavily on the relationship between a\, . . . , a m , b±, . . . , b n . 

The fact that the counting function T(a, b) is of a somewhat more delicate nature naturally has 
computational consequences. We feel that providing any general results on this function would go 
beyond the scope of this article and will hopefully be the subject of a future project. On the other 
hand, we adjusted our algorithm to compute values of T(a, b) for three (fixed) pairs (a, b), which 
have been previously computed by Mount [17] and DeLoera and Sturmfels [9]. The first example is 

T((3046, 5173, 6116, 10928), (182, 778, 3635, 9558, 11110)) = 
23196436596128897574829611531938753 

The authors reported their computation took 20 minutes (Mount) or 10 minutes (DeLoera-Sturmfels) 
[9, 17]. We computed this number in 0.2 seconds, based on Theorem 5. A similar phenomenon 
happens with 

T((338106, 574203, 678876, 1213008), (20202, 142746, 410755, 1007773, 1222717)) = 
316052820930116909459822049052149787748004963058022997262397 

and 

T((30201, 59791, 70017, 41731, 58270), (81016, 68993, 47000, 43001, 20000)) = 

24640538268151981086397018033422264050757251133401758112509495633028 

which reportedly took 35 minutes/10 days with the DeLoera-Sturmfels algorithm [9], 0.3/2.9 sec- 
onds with ours. These timing comparisons ignore differences in machine speeds and implementation 
of the algorithms, but suggest that our methods are considerably more efficient. 



Acknowledgements. We would like to thank the referees for carefully reading through our paper 
and software, and for their many helpful comments. 
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